# Import data from .RData
#load("data/PSZ_cont_snotel.RData")
load("data/cont_snotel.RData")
# Count the number of SNOTEL Sites per Ecoregion
eco_sntl_sites <- cont_snotel %>%
group_by(US_L3NAME) %>%
summarize(eco_sntl_sites = n_distinct(site_id)) %>%
filter(eco_sntl_sites >= 5)
# All SNOTEL Sites, Calculate spring days with new SWE and new SWE depth
# filter by date (March 1st - 60/61, March 15th - 74/75, April 1st - 91/92, April 15th - 105/106, August 1st - 243/244)
# Filter data by date from 2000 to present and from March 1st to August 1st
# group by year and site
start_year = 2000
# start_day = 60 #March 1
# cont_snotel_spring_mar1 <- cont_snotel %>%
# mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
# mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
# filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243,
# yday >= start_day+1 & yday <= 243+1)) %>%
# select(-c("X", "network", "description", "start", "end"))
#
# start_day = 74 #March 15
# cont_snotel_spring_mar15 <- cont_snotel %>%
# mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
# mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
# filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243,
# yday >= start_day+1 & yday <= 243+1)) %>%
# select(-c("X", "network", "description", "start", "end"))
start_day = 91 #April 1
cont_snotel_spring_apr1 <- cont_snotel %>%
mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent),
snow_water_equivalent-lag(snow_water_equivalent), 0),
newSWE = ifelse(temperature_min > 0, 0, newSWE)) %>%
filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243,
yday >= start_day+1 & yday <= 243+1)) %>%
select(-c("X", "network", "description", "start", "end"))
# start_day = 105 #April 15
# cont_snotel_spring_apr15 <- cont_snotel %>%
# mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
# mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent)+3, snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
# filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243,
# yday >= start_day+1 & yday <= 243+1)) %>%
# select(-c("X", "network", "description", "start", "end"))
# # Determine number of days with increased SWE per spring at each site (YEARLY - March 1)
# cont_snotel_site_yr_mar1 <- cont_snotel_spring_mar1 %>%
# group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
# summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>%
# merge(., eco_sntl_sites, by = c("US_L3NAME")) %>%
# mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "3/1")
#
# # Determine number of days with increased SWE per spring at each site (YEARLY - March 15)
# cont_snotel_site_yr_mar15 <- cont_snotel_spring_mar15 %>%
# group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
# summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>%
# merge(., eco_sntl_sites, by = c("US_L3NAME")) %>%
# mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "3/15")
# Determine number of days with increased SWE per spring at each site (YEARLY - April 1)
cont_snotel_site_yr_apr1 <- cont_snotel_spring_apr1 %>%
group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>%
merge(., eco_sntl_sites, by = c("US_L3NAME")) %>%
mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "4/1")
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude', 'elev', 'US_L3CODE', 'US_L3NAME'. You can override
## using the `.groups` argument.
# # Determine number of days with increased SWE per spring at each site (YEARLY - April 15)
# cont_snotel_site_yr_apr15 <- cont_snotel_spring_apr15 %>%
# group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
# summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>%
# merge(., eco_sntl_sites, by = c("US_L3NAME")) %>%
# mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "4/15")
# Merge the DF's into one large datafame for plotting
# cont_snotel_site_yr <- rbind(cont_snotel_site_yr_mar1, cont_snotel_site_yr_mar15, cont_snotel_site_yr_apr1, cont_snotel_site_yr_apr15)
# Determine average number of days with increased SWE per spring at each site in the PSZ (SITE - April 1)
cont_snotel_site_apr1 <- cont_snotel_site_yr_apr1 %>%
group_by(site_id, site_name, state, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY, eco_sntl_sites, eco_sntl) %>%
summarise(mean_spring_days = mean(newSWE_days), med_spring_days = median(newSWE_days),
mean_newSWE = mean(newSWE), med_newSWE = median(newSWE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude', 'elev', 'US_L3CODE', 'US_L3NAME', 'L3_KEY',
## 'eco_sntl_sites'. You can override using the `.groups` argument.
# # View data
# mapview(PSZ_cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "med_spring_days",
# layer.name = "AVG Sping Days (PSZ Sites)" , crs = 4326, grid = FALSE) +
# mapview(eco_L3_4326, zcol = "L3_KEY")
# Spring SWE increases (After April 1)
spring_snow_days <- ggplot(cont_snotel_site_yr_apr1, aes(x = reorder(US_L3CODE,newSWE_days,na.rm = TRUE), y= newSWE_days,
fill = reorder(eco_sntl,newSWE_days,na.rm = TRUE))) +
geom_boxplot() +
labs(x = "EPA Level III Ecoregion", y = "Days of Increased SWE", fill = "") +
scale_y_continuous(breaks = seq(0,50,10), limits = c(0,50), expand = c(0.01, 0.01)) +
scale_x_discrete(labels = function(eco_sntl) str_wrap(eco_sntl, width = 28)) +
theme_bw()
spring_snow_days
ggplotly(spring_snow_days)
#ggplot_build(spring_snow_days)$data
# Pairwise comparison of all sites
stat.test <- cont_snotel_site_yr_apr1 %>%
mutate(US_L3CODE = as.numeric(US_L3CODE)) %>%
arrange(US_L3CODE) %>%
pairwise_t_test(newSWE_days ~ US_L3CODE, p.adjust.method = "bonferroni") %>%
add_significance(cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns"))
level1 = c(4,5,9,11,13,15,16,17,18,19,21,23,41,77,80)
level2 = c(5,9,11,13,15,16,17,18,19,21,23,41,77,80)
pairwise_plot <- ggplot(stat.test, aes(x = factor(group1, levels = level1),
y = factor(group2, levels = level2),
fill = p.adj.signif))+
geom_tile(col = "black") +
scale_fill_brewer(labels = c("1e-04", "0.001", "0.01", "0.05", "ns"),
palette = "Reds", direction = -1)+
labs(y = "EPA Level III Ecoregion", x = "EPA Level III Ecoregion", fill = "Sig. Level") +
theme_bw()
pairwise_plot
ggplotly(pairwise_plot)
#test <- ggplot_build(pairwise_plot)$data[[1]]
# boxplot and pairwise
plots <- plot_grid(
spring_snow_days + theme(legend.position="none"),
pairwise_plot,
align = 'vh',
labels = c("a.", "b."),
hjust = -0.4,
nrow = 1)
## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.
plots
legend <- get_legend(spring_snow_days + guides(fill = guide_legend(nrow = 5)) +
theme(legend.direction = "horizontal",
legend.justification="center",
legend.box.just = "bottom"))
plot_grid(plots, legend, ncol = 1, rel_heights = c(1, .3))
# Elevation vs new SWE count (March 1)
elev_newsnow <- ggplot(cont_snotel_site_apr1, aes(x = med_spring_days, y= elev,
color = reorder(eco_sntl,med_spring_days,na.rm = TRUE))) +
geom_point() +
geom_smooth(method = "lm")+
stat_smooth(method = "lm", col = "black")+
labs(x = "Median Days of Increased SWE", y = "Elevation (m)", color = "") +
scale_x_continuous(breaks = seq(0,30,10), limits = c(0,30), expand = c(0.01, 0.01)) +
scale_y_continuous(breaks = seq(500,4000,1000), limits = c(500,3750), expand = c(0.01, 0.01)) +
scale_color_manual(values=c("#F8766D", "#C99800", "#E58700", "#A3A500", "#6BB100",
"#00BA38", "#00BF7D", "#00C0AF", "#00B0F6", "#619CFF",
"#B983FF", "#00BCD8", "#E76BF3", "#FF67A4", "#FD61D1")) +
theme_bw()
elev_newsnow
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggplotly(elev_newsnow)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
#ggplot_build(elev_newsnow)$data
plot <- plot_grid(
spring_snow_days + theme(legend.position="none"),
elev_newsnow + theme(legend.position="none"),
align = 'vh',
labels = c("a.", "b."),
hjust = -1,
nrow = 1)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
plot
legend <- get_legend(spring_snow_days + guides(fill = guide_legend(nrow = 5)) +
theme(legend.direction = "horizontal",
legend.justification="center",
legend.box.just = "bottom"))
plot_grid(plot, legend, ncol = 1, rel_heights = c(1, .3))
#Determine the number of SWE events above a temperature threshold
summary_stats <- cont_snotel_spring_apr1 %>%
#filter(site_id == "1000" & (date >= "2001-04-01" & date <= "2001-04-15")) %>%
summarise(n_above_mean = sum(newSWE > 0 & temperature_mean > 0, na.rm = TRUE),
n_above_min = sum(newSWE > 0 & temperature_min > 0, na.rm = TRUE),
n_below_mean = sum(newSWE > 0 & temperature_mean <= 0, na.rm = TRUE),
n_below_min = sum(newSWE > 0 & temperature_min <=0, na.rm = TRUE),
n_newSWE = sum(newSWE > 0, na.rm = TRUE),
n = n())
above_min_pct = summary_stats$n_above_min/(summary_stats$n_newSWE)